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Abstract 



We present a numerically tractable method to solve exactly the evolution of a N boson system with 
binary interactions. The density operator of the system p is obtained as the stochastic average of particular 
operators |^'i)(4 , 2| of the system. The states ^1,2) are either Fock states \N : $1,2) or coherent states 
|coh : 0i,2) with each particle in the state (j>i 2(x). We determine the conditions on the evolution of 
4>i,2 -which involves a stochastic element- under which we recover the exact evolution of p. We discuss 
various possible implementations of these conditions. The well known positive P-representation arises as 
"-^J ] a particular case of the coherent state ansatz. We treat numerically two examples: a two-mode system 

and a one-dimensional harmonically confined gas. These examples, together with an analytical estimate 
of the noise, show that the Fock state ansatz is the most promising one in terms of precision and stability 
of the numerical solution. 



Pacs: 03.75.Fi, 05.10.Gg, 42.50.-p 

> 

m 1 Introduction 

m ■ 
o 

0> Since the experimental realization of the first atomic gaseous Bose-Einstein condensates a few years ago 
[||, H], ||, |4j, the physics of dilute Bose gases has been considered with a renewed interest. One fascinating 
£3 ■ aspect of these new systems is the possibility to accumulate in a single quantum state a large fraction of the 
atoms confined in a trap. 



At very low temperature, a simple theoretical description of the dynamics of these systems is obtained by 
C . neglecting the uncondensed atoms, and by considering the wave function of the condensate, which obeys a 
Schrodinger equation with a non linear term originating from the mean-field interactions between the atoms. 
Such an approach neglects two- and more-particle correlations and is valid under a weak-interaction condition 
which is usually stated in terms of the density n and the scattering length a of the gas as (na 3 ) 1 / 2 <C 1. 
Current gaseous condensates satisfy such a condition. Nevertheless effects beyond the Gross-Pitaevskii 
' equation may be considered at zero temperature; also finite temperature phenomena are not accounted for 
by the pure state mean field approach. 

More complex theories have been developed in order to cope with effects beyond the Gross-Pitaevskii 
equation: Bogoliubov's approach takes into account the next term in the (na 3 ) 1//2 expansion |5|, Also 
quantum kinetic theories have been developed to study the formation of the condensate and to include the 
effect of the non-condensed particles J?], ||, Unfortunately the corresponding calculations are quite heavy 
for 3D non-homogeneous systems such as trapped gases, and this constitutes a first limitation to the use of 
these methods. Also approximations used in some of these mean field theories are not under rigorous control, 
making it difficult to assess their domain of validity (for a review see e.g. |l(J). Therefore a computational 
scheme capable to provide exact results can have a great importance both from a purely theoretical point 
of view and for a quantitative analysis of experimental data. 
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When the Bose gas is at thermal equilibrium such an exact numerical calculation of the properties of the 
gas is available, using the Quantum Monte-Carlo techniques, based on Feynman path integral formulation 
of quantum mechanics fTl| , |l2]| . The aim of this paper is to present an alternative exact and numerically 
tractable solution to the problem of the interacting Bose gas, a method not restricted to the case of thermal 
equilibrium but which allows for the study of the dynamics of the gas. The method is based on a stochastic 
evolution of Hartree states, in which all atoms have the same wave function, these Hartree states being 
either Fock states (fixed number of atoms) or coherent states. As a particular case of this solution with 
coherent states, we recover the stochastic scheme corresponding to the evolution of the density operator 
of the system in the positive P-representation [jl3| , |l4| , |i~5| , [l6}| . This evolution is known to lead to strong 
unstability problems, which fortunately do not show up for other implementations of the present method. 

The outline of the paper is the following: In section ^, we present the stochastic formulation of the 
evolution of these Hartree states which, after average over the stochastic component, leads to the exact 
evolution. Section [3| is devoted to the presentation of two particular schemes implementing this stochastic 
formulation. We first present a 'simple' scheme, which minimizes the statistical spread of the calculated 
iV-atom density matrix. We also investigate a more elaborate scheme in which the trace of the calculated 
density matrix remains strictly constant in the evolution. With this constraint, we recover for coherent states 
the known stochastic simulation associated with the positive P-representation Finally we investigate 
in sections |] and || two examples, a two-mode model system and a one-dimensional Bose gas respectively. 
These examples illustrate the accuracy and the limitations of the method. Generally speaking we find that 
the 'simple' scheme simulations are only limited by the computation power: the number of realizations 
needed for a good statistical accuracy increases exponentially with time for the simulation with Fock states. 
On the contrary the simulations with constant trace are subject to divergences of the norms of the stochastic 
wave functions in finite time, a phenomenon already known for coherent states in the context of the positive 



P-representation [17]. 



2 Stochastic formulation of the TV-boson problem using Hartree func- 
tions 

2.1 Model considered in this paper 

The Hamiltonian of the trapped interacting Bose gas under exam can be written in terms of the Bose field 
operator i&(x) as: 

H = Jdx&(x)h $(x) + ~ J J dxdx'i>\x)¥(x')V(x-x')i>(x')i>{x) (1) 

where x is the set of spatial coordinates of a particle, ho = — +V cyi t{x) is the single particle Hamiltonian 
in the external confining potential V ex t and where interactions are assumed to occur via a two-body potential 
V(x-x'). 

In practice we consider the dilute gas and the low temperature regimes, which correspond respectively 
to n\a\ 3 <C 1 and \a\ <C A for a three-dimensional problem (A = /i/(27rm/cBT) 1//2 is the thermal de Broglie 
wavelength). The true interaction potential can then be replaced by a simpler model potential leading to 
the same scattering length a provided that the range b of this model potential is much smaller than the 
healing length £ = (87rna) -1 / 2 and than A (5|, |(|. This ensures that the physical results do not depend on b. 
For simplicity we will use here repulsive Gaussian potentials corresponding to a positive scattering length 
a > 0. 
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2.2 A stochastic Hartree Ansatz with Fock states 



From a mathematical point of view, the exact evolution of the iV-body density matrix p can be obtained 
from the Hamiltonian (|l]) using the quantum-mechanical equation of motion 

p(t) = ^[H,p{t)] (2) 

but any concrete calculation is impracticable even for moderate particle numbers N, due to the multi-mode 
nature of the problem leading to a huge dimensionality of the N— body Hilbert space. 

For this reason approximate theories have been developed in order to get useful results at least in some 
specific ranges of parameters; the simplest one is the so-called mean-field theory, in which the iV-particle 
density matrix is approximated by a Fock state Hartree ansatz 

p(t) = \N:</>(t))(N:<j>(t)\. (3) 

The evolution of the normalized condensate wave function (j) is determined using either a factorization 
approximation in the evolution equation for the field operator || or a variational procedure IS]. The result 
is the well-known mean-field equation 



in ^L = (-^- + V ext (x)) <j>(x) + (N - 1) (J dx' V(x - x') \<j>(x')\ 

For an interaction potential V(x — x') modeled by a contact term g5(x — x') (where g = 4irfi 2 a/m in a three- 
dimensional problem) it reduces to the Gross-Pitaevskii equation commonly used to analyze the dynamics 
of pure Bose-Einstein condensed gases. 

A first attempt to improve the accuracy of the Hartree ansatz ([|) is to allow for a stochastic contribution 
dB in the evolution of the macroscopic wave function (f>: 

(f,{t + dt) = <f>(t) +F dt + dB . (5) 

In all this paper the noise dB is treated in the standard Ito formalism IS]: it is assumed to have a zero mean 
dB = and to have a variance dB 2 oc dt; a deterministic contribution is given by the "force" term Fdt. In 
this framework, the iV-body density matrix would result from the stochastic mean over noise or, in other 
terms, from a mean over the probability distribution V{4>) in the functional space of the wave functions (j>: 

p(t)l(\N:<f>(t)){N:<f>{t)\) = fv<t>V(<t>) \N : <t>(t))(N : cf>(t)\. (6) 

\ / stoch J 

An immediate advantage of this prescription over the pure state ansatz Eq.(0) is that it could deal with 



finite temperature problems [20|. However as shown in ^2^, the simple generalization Eq.(|5|) of the Gross- 
Pitaevskii equation cannot lead to an exact solution of the iV-body problem [^TJ. Therefore we have to 
enlarge the family of dyadics over which we expand the density operator; more precisely we use Hartree 
dyadics in which the wave functions in the bra and in the ket are different: 

a(t) = \N :Mt))(N :Mt)\ ■ (7) 

The two wave functions 4>\{x) and 4>2(x) are assumed to evolve according to Ito stochastic differential 
equations: 

4> a (t + dt) = 4> a (t) + F a dt + dB a (a = 1,2). (8) 
The expansion Eq.(^|) is then replaced by 

p(t) = (\N : Mt))(N : Mt)\) stQch = J jvfcVfaV^fa) \N : 0i(t) ){N : 2 (i)|. (9) 
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We will see in the following that within this extended Hartree ansatz one can find a stochastic evolution for 
01,2 reproducing the exact time evolution. 

Actual calculations (see §|4] and §||) will be performed with a Monte-Carlo technique, in which the 
evolution of the probability distribution V is simulated by a large but finite number M of independent 
realizations 4>{\{t), i = 1, . . . ,J\f. At any time the (approximate) density matrix p is given by the mean 
over such an ensemble of wave functions: 

1 N 

p(t)^^\N:4\t))(N:4\t)\. (10) 
1=1 

The expectation value of any operator O is thus expressed by: 

1 N 

[ 6 )^JfT,( N: ^\t)\0\N : 4\t)) . (11) 

i=l 

For an Hermitian operator one can equivalently consider only the real part of this expression since the 
imaginary part is vanishingly small in the large M limit. 

Consider as an example the one-particle density matrix of the gas, usually defined as: 

p (1) (x,x') = (&(x')$(x)). (12) 



Inserting in this expression our form of the complete density matrix (10), we obtain the simple result 

pV{x,x') = N (Mx)tiW (^i) 1 *- 1 ) u (13) 

\ / stoch 

from which it is easy to obtain the spatial density n(x) = p^(x, x) and the correlation function g^\x, x') = 
p^\x, x') / {n{x)n{x')) 1 / 2 . Also, the condensate fraction can be obtained from the largest eigenvalue of 
pW(x,x'). 

Remarks: 

1. The desired stochastic evolution, which has to satisfy Tr[p] = 1, cannot preserve the normalization of 
01,2 to unity; we can write indeed 

Tv[p(t)] = ((Mt)\Mt)) N ) „ = 1 (14) 

\ / stoch 

which for \cpi) ^ \cj>2) imposes \\4>i\\ ||02|| > 1- 

2. The expansion Eq.(Q) is always possible. Using the identity 

1 M 

Id N = lim — \N : tp^}(N : ipV)\ (15) 

where the functions ij)V> have a uniform distribution over the unit sphere in the functional space, we 
obtain 



p= lim V \N :^ jl) )(N :^ j2) \(N :^ jl) \p\N :^ j2) ). (16) 

n j2=i 

We write the matrix elements (N : ip^^\p\N : ip^) as C(^j 2 ) an d we se t <Px'^ = V'^CO'iJa) an£ l 
_ ^0' 2 ) £* Putting N = A4 2 and reindexing (ji,j2) as a single index i we recover the 
expansion Eq.(|l0|). Note that this expansion is not unique and does not have the pretension to be the 
most efficient one. For instance if the system is initially in a Hartree state \N : 0o)> such a procedure 
is clearly not needed since one has just to set (f>\ (t = 0) = = 0) = 4>o- This will be the case of 

the numerical examples in sections Q and [El. 
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2.3 Stochastic evolution of a Fock state Hartree dyadic 



In this subsection we calculate the stochastic time evolution during an infinitesimal time interval dt of the 
dyadic a(t) given in Eq. (]?]). This will be used later in a comparison with the exact master equation. 

After dt, the dyadic a has evolved into: 

a(t + dt) = \N : fa + <% )(N : 4>i + dfa\, (17) 

where d<f>\ and dfa, defined according to (||), contain both the deterministic contribution F a dt and the 
stochastic one dB a . Splitting each contribution into a longitudinal and an orthogonal component with 
respect to <f> a and isolating a Gross-Pitaevskii term in the deterministic contribution, we can write: 



dB a {x) = 4> a {x) d~f a + dB^(x) 
F a (x) = F^ p (x) + \ a <j) a {x) + F^(x). 

Our choice of the Gross-Pitaevskii term is the following one: 



(18) 
(19) 



F? p (x) 



1 

ih 

1 

ih 



h + 



^ j dx V(x - x) \ (j) a {x')\ 
WaW J 

(N- 1) (<f> a <f>a\V\<l>a<l>a, 

<p a \X). 



4>a(x) 



(20) 



The first term gives the standard Gross-Pitaevskii evolution, including the kinetic term, the potential energy 
of the trap and the mean-field interaction energy; the second term, which arises naturally because we are 
considering Fock states (rather than coherent states as commonly done) takes into account the difference 



between the total mean-field energy per particle of the condensate and its chemical potential [i [22] 



We split the field operator in its longitudinal and transverse components, keeping in mind that the wave 
functions <b n are not of unit norm: 



W{x) 



,t 



with 



The relevant bosonic commutation relations then read: 

= Haf and [a^ , (ar)] = 0. 
We will also need the projector Q a onto the subspace orthogonal to 4> a : 



Q^mx,x',...)] = ^(x,x',...) 



dy 4>* a {y) ip(y,x', . 



(21) 

(22) 
(23) 
(24) 



This projector arises in the calculation as we have introduced a component of the field operator orthogonal 
to 4> a . Using J dx <f> a (x) 5^fa(x) = we shall transform integrals involving 8^S>a(x) as follows: 



dxi>(x,x\...)8$l L (x) = dx Q^[t/j(x,x' , . . .^^(x). 



(25) 



Inserting these definitions in (|l~7|) the expression for a at time t + dt can be written as 
a(t + dt) - a(t) = S[ 0) \N :</>x){N: <fo\ +e.c. 

+ j dx S[ 1] (x)5$\(x)\N - 1 : fa){N : <fo\ + e.c. 

+ / [dxdx's{ 2) (x,x f )5y{(x)5y{(x')\N-2:fa)(N:<f)2\+e.c. 

+ I fdxdx' S^(x,x')S^\(x)\N - 1 : fa)(N- 1 : <h\^2{x') 



(26) 
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where the notation e.c. stands for the exchanged and conjugate of a quantity, that is the complex conjugate 

(i) 

of the same quantity after having exchanged the indices 1 and 2. The explicit expressions for the Sa are: 



Si 



(0) 



N 



dt + NXi dt + Nd'j! 



N(N - 1) 



N 2 



2 dli + ~Y d ^ d 72 



(27) 



S[ 1) ( x ) = VN {qS x) [F? p (x)] dt + F^(x)dt + dB^s) + (iV - l)d 7 i dB^x) + iV dB^(x) d 72 *}(28) 



S ( ^(x,x') = vESLJldB^x) dBj-(x') 
5 (1 ' 1} (x,x') = N dBi(x) dBt{x'). 



(29) 
(30) 



Analogous expressions for S% , are obtained by exchanging the indices 1 and 2. In the next 

subsection, we evaluate the exact evolution of the same dyadic during a time interval dt, so that we can 
determine the constraints on the force and noise terms entering into these equations. 



2.4 Exact evolution of a Fock state Hartree dyadic 

To make the stochastic scheme described in the previous sections equivalent to the exact dynamics as it is 
given by ([[]), the final result of the previous subsection (|26|)-(j30l) has to be compared with the exact evolution 
of the density matrix a(t). Consider a dyadic a = \N : 0i )( N : 02 1 at time t; according to the equation of 
motion (^), after an infinitesimal time step dt it has evolved into: 

ait + dt) = ait) + -4 CHait) - a(t)H) 
in 

= a(t) +E{ 0) \N : fa){N : <fo\ + e.c. 
+ fdx e[ 1 \x)S^\{x)\N - 1 : 0i)(JV : 2 | +e.c. 

-!- / I dx dx' (x,x')5^\(x)5^\(x')\N - 2 : cj) 1 )(N :0 2 |+e.c. 



(i) 

where the E& are given by 

,(o) _ Ndt 



E 



ih 



(0i|fc o |0i) , (at- l) (^i^il^i^i^! 



+ 



EP(x) 



dWN {x) 
in 



2 

(N-l) 



dx' V(x — x') |0i(x / )| 2 ) 4>i(x) 



E ?\ x , x >) = 1} Q[ x) Q^[V(x-x')Mx)Mx% 

Analogous expressions for E^\ E^\ are obtained by exchanging the indices 1 and 2. 



(31) 

(32) 
(33) 
(34) 



2.5 Validity conditions for the stochastic Fock state Hartree ansatz 



The similarity of the structures of ( |26|) and (31) suggests the possibility of a stochastic scheme equivalent 



to the exact evolution: to achieve this, it is necessary to find out specific forms of deterministic (19) and 
stochastic ( |l~8| ) terms for which the mean values of the equal the E& : 



s£\x)=E$\x) 



Set 3?') — -^a ^ ^ ) 

S^ix,^) = 0. 



(35) 

(36) 

(37) 
(38) 
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In terms of the different components, these conditions can be rewritten as: 

(N-l) 



+ Ndjidj^ = 



From the last equation (38), it follows immediately why independent bras and kets are needed in the 
ansatz ([/J): in the case <f>i = (f>2 = <fi such a condition would in fact lead to a vanishing orthogonal noise and 
finally to the impossibility of satisfying (| 



(39) 

(40) 
(41) 

(42) 
(43) 



(Ai + X* 2 ) dt + 



F{-(x) dt+(N- l)dB^[x) d 7 i + NdBj-(x) d^ = 
F 2 L (x) dt+(N- l)dB 2 L (x)dj 2 + NdB^(x)djl = 
QaQ X a [V{x - x')(j) a (x)(j) a (x')] 



dB^x)dB^(x') 



dt 
Hi 



dB^{x)dB 2 L *{x') = 0. 

As we shall discuss in detail in §||, several different stochastic schemes can be found satisfying (pii|)-(|4"3|); 
each of them gives an evolution identical in average to the exact one, but the statistical properties can be 
very different. 



2.6 A stochastic Hartree ansatz with coherent states 

Up to now we have worked out the case of a Fock state ansatz \N : (j>i)(N : (f> 2 |. Actually coherent states 
rather than Fock states are generally used, both in quantum optics and in condensed matter physics. We 
now show that our stochastic procedure also applies with a coherent state ansatz of the form: 

a(t) = II(t) | coh : 0i )( coh : <£ 2 |, (44) 

with 

|coh:«=exp(^y (faWx)# . W ) |vac) (45) 

where N is the mean number of particles. We have included here a prefactor II(i) which was absent in the 
case of the Fock state ansatz Eq. (0) ; in the Fock state case indeed such a prefactor could be reincluded into 
the definition of 4>\ and cp2- The wave functions (j>a{x) and the prefactor factor II evolve according to Ito 
stochastic differential equations 



Splitting the field operator as 
and using 



d(j) a = F a dt + dB a 

dU = fdt + db. (46) 

^(x) = N 1/2 <j) a {x) + 5$ a (x) (47) 

8$ a (x) | coh :<t> a } = (48) 



we find that the equivalence of the stochastic scheme and the exact evolution translates into the set of 
conditions: 

/ = (49) 

Fi{x)dt + l-dbdBAx) = ^hofatx) (50) 
II in 

F 2 (x)dt + ^-db* dB 2 {x) = ^rhoMx) (51) 
11* in 

dt. 



dBJx)dBJx') = —V(x - x')4> a (x)cl> a (x') (52) 
in 

dB 1 {x)dB^(x l ) = 0. (53) 
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As we shall see in §|3[ such conditions are satisfied by several stochastic schemes. Very remarkably, the 
stochastic evolution deduced from the positive P-representation [23] arises naturally as one of them. 

Within this coherent state ansatz the one-particle density matrix is evaluated using 

pW(x,x') = NU 1 {x)4>l{x')Ii{t) exp(N {falfa))) . (54) 

\ / stoch 

In a practical implementation of the simulation it turns out to be numerically more efficient to represent 
n(t) as the exponential of some quantity 

n(t) = e Ns ^ (55) 
and to evolve S(t) according to the stochastic equation 

dS = - J dx [dB!(x)<l>*(x) + dB%(x)fa(x)] 

A,// f dx [ dx'Vix-x^llM^lMxW-lM^flMx')] 12 ]- (56) 



2ih 

3 Particular implementations of the stochastic approach 

In the previous section we have derived the conditions that a stochastic scheme has to satisfy in order to 
recover the exact evolution given by the Hamiltonian (|l|); in the case of the Fock state ansatz (0), we get to 
the system (|39"1)-([43|), while in the case of the coherent state ansatz ( fl4| ) we get to the conditions (^g|)-([53]). 
As the number of these equations is actually smaller than the number of unknown functions there is by 
no mean uniqueness of the solutions, that is of the simulation schemes. We need a strategy to identify 
interesting solutions. 

We therefore start this section by considering various indicators of the statistical error of the simulation 
(§|3,l|) which can be used as guidelines in the search for simulation schemes. These indicators are defined 
as variances of relevant quantities which are conserved in the exact evolution but which may fluctuate in 
the simulation. We show that these indicators are non decreasing functions of time; attempts to minimize 
the time derivative of a specific indicator will lead to particular implementations of the general stochastic 
method, such as the simple scheme (§|3.2|) and the constant trace scheme (§|3.3|). 



3.1 Growth of the statistical errors 

The first indicator that we consider measures the squared deviation of the stochastic operator a it) from the 
exact density operator p{t): 



A(t) = h?[{o\t) - p{t)){o{t) - p{t))\ 

\ / stoch 

W(i)a(t)A -Tr[p(i) 2 ]. (57) 



stoch 



We now show that A(t) is a non-decreasing function of time. When the stochastic scheme satisfies the 
validity conditions derived in the previous section, we can write the stochastic equation for a as: 

da = \n,cr} + da s (58) 
in 

where da s is a zero- mean noise term linear in dB a (and db for the coherent state simulation). In the case of 
simulation with Fock states it is given by 

da s = N 1/2 ^J dx dBi(x)¥ \x)\N - 1 : <j>i)(N : <fa\ + J dx dBl(x)\N : 0i)(JV - 1 : </> 2 |*(a;)| ■ (59) 
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In the case of simulation with coherent states it is given by 



da s = db\coh : </>i)(coh : 02 \ 

+ iV 1/2 n<! I dxdB^x)^ (x)\coh :0i)(coh : fa\ + I dx dB^(x)\coh :0i)(coh : ^a|*(x) .(60) 



We calculate the variation of A during dt, replacing cr by cr + da in Eq.([57]) and keeping terms up to 
order dt. Using the invariance of the trace in a cyclic permutation and averaging over the noise between t 
and t + dt we finally obtain 



dA 



Tr[da\da s 



stoch 



(61) 



which is a positive quantity. Minimization of this quantity is the subject of §0. Physically dA > means 
that the impurity of the stochastic density operator a always increases in average, while the exact density 
operator has a constant purity Tr[p 2 ]. 

The second kind of indicator that we consider measures the statistical error on constants of motion of the 
exact evolution. Consider a time independent operator X commuting with the Hamiltonian. The stochastic 
evolution leads to an error on the expectation value of X with a variance given by the ensemble average of 



Ax(t) 



Ti[Xa(t)] - Tr[Xp(t)] 



stoch 



Tv[Xa(t)} 



stoch 



Tv[Xp(t)} 



(62) 



From Eq.(p8|) we obtain the variation after a time step dt of the expectation value of X along a stochastic 
trajectory: 



d(Tr[Xa]) = -^Tr (X[H, o"]) + Tr[Xda s ]. 
in 



(63) 



Using the invariance of the trace under cyclic permutation and the commutation relation [7i, X] = we find 
that the first term in the right hand side of Eq. vanishes so that 



dA x 



Tr [Xda s 



2 ) 

I stoch 



(64) 



a quantity which is always non-negative. 



Using expression (|64j) one can 'design' simulations preserving exactly the conserved quantity, the con- 
straint to meet being Tr[Xdo" s ] = 0: for instance in the Fock state simulation, one first chooses the transverse 
noises dB^ satisfying Eqs,( f42| , f4li| ); then one simply has to take for the longitudinal noise of <f>i: 



7^{( N ■ <P2\X\N : ^i))- 1 J dx dBj-(x)(N : (P 2 \X5^\(x)\N - 1 : fa) 



(65) 



and a similar expression for d"f2', finally the force terms F Q are adjusted in order to satisfy Eq. (39-41). As 
natural examples of conserved quantities one can choose X = 1 or X = TC; the former case is discussed in 
detail in SpT 



3.2 The simple schemes 

These schemes are characterized by the minimization of the incremental variation of the statistic spread 
of the A^-particle density matrix a(t), a spread that we have already quantified in Eq.([57]) by A(t). To be 
more specific we assume that we have evolved a dyadic up to time t, and we look for the noise terms that 
minimize the increase of Trfa^o"] between t and t + dt. 
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3.2.1 Simulation with Fock states 



In the case of the Fock state ansatz, we calculate explicitly the variation of Tr^er] from Eq.(]59|) and we 
get: 



iVTr[crt(T] 



N\d ll +d 1 *\* + 



a=l,2 



dx \dB^{x)\ 2 + [dr/i + d^ + ex. 



(66) 



We now look for the noise terms d^y a and dB^ minimizing this quantity under the constraints Eqs.([3S|-43). 
We first note that we can choose d"fi = d^2 = without affecting the transverse noises, as shown by 



Eqs.(3£-43): the correlation function of the transverse noises do not involve the dj a , and we can accom- 
modate for any choice of d^ a by defining appropriately the force terms , \ a • In the particular case 
defining our simple scheme we take all these force terms equal to zero. Note that the choice of vanishing 
dj's immediately leads to a vanishing noise term in Eq.(|66|). 



Secondly the terms involving the transverse noise in Eq.(^) are bounded from below: As the modulus 
of a mean is less than the mean of the modulus, we have 



dB±(x)dB£{x) < \dB±(x)\ 2 



with the left hand side of this inequality fully determined by condition Eq. 

We have found for the transverse noise a choice which fulfills Eqs.(|4^ 
inequality Eq.(|67|): 



dB^{x) 



1/2 



nO) 



dk 



(27r) rf / 2 



V(k) 



1/2 



(67) 



and which saturates the 



e ikx e id a (k) 



where d is the dimension of position space, V{k) is the Fourier transform of the model interaction potential, 
supposed here to be positive: 



V(k) = / dx V(x)e 



-ikx 



(69) 



The phases 6 a have the following statistical property: 



je a (k) e i8 a (k>) _ fiffc _|_ y\ 



(70) 



and Ox, 02 are uncorrelated. In practice for half of the /c-space (e.g. k z > 0) 6 a (k) is randomly chosen between 
and 2tt; for the remaining fe's we take a (—k) = —6 a (k). One can then check that this choice for the 
transverse noise reproduces the correlation function Eqs.(|4^j43]). 



We show now that the implementation ( p8|) saturates the inequality Eq. (p7|) , so that it leads to the 
minimal possible value for dTrlcro"] within the validity constraints of the stochastic approach. We calculate 



explicitly the right hand side of Eq.(67) 



dt 



V(0) 



r ltd (V )\\rf V . 

dy n , M2 V(x -y) + 

\\<Pa\\ 



\v\ 



(71) 



where Q* projects onto the subspace orthogonal to 0* and where we have used the positivity of the Fourier 
transform V of the model interaction potential. The left hand side of Eq.(pTj) is calculated using Eq.(f42|): 



dt 



dB^xf = -ct>i{x) 



V(0) 



2 / d y m M2 v ( x -y) + - 



(72) 
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As the expressions between square brackets in Eqs.( fF^|7l| ) are real positive we deduce the equality in Eq.(^ 
We can now calculate explicitly the variation of Trfa^cr] by integrating Eq.(|7l|) over x: 



cZTr[o-t(j] _ dt 
iVTrfcrto-] ~ 7T 



2V(0) - ' 

a=l,2 



W\ 



1100 



(73) 



This expression is particularly useful since it allows one to derive an upper bound on the increase of Trfa^cr]: 
as we assume here a positive Gaussian model potential V(x — x') the matrix element {(fiai^a \V\4> a -,4>a) is 
positive so that the right hand side of Eq.([73]) is smaller than 2V(0)dt/h. After time integration, using 
Eq.(|57|) and the fact that the trace of the squared density operator p 2 is a constant under Hamiltonian 
evolution we can deduce an upper bound on the squared statistical error A(i): 



A(i) + Tr[p 2 ]< [A(0) + Tr[p 2 ]]e 



2NV(0)t/h 



(74) 



Note that it involves the model dependent quantity V(0) and not only the physical parameters of the problem 
such as the chemical potential or the scattering length. It may be therefore important to adjust the model 
interaction potential V{x — x') in order to minimize the growth of the statistical error for given physical 
parameters. 

To summarize the proposed simple scheme has several noticeable properties. The deterministic force 
acting on the </> Q 's is simply the mean field contribution, so that the whole correction to the mean field 
evolution is provided by the transverse noises dB^ . Also the evolutions of the two states <f> a are totally 
independent from each other. 



3.2.2 Simulation with coherent states 

In the case of the coherent state ansatz, an explicit calculation of dTr[(T^cj] from Eq.(^) gives: 
dTrftrtff] 



Tr[crta] 



db 



+ N / dx dBi(x)(f>l{x) + dB*(x)(f)2(x] 



+ N Y I dx \dB 
o=l,2 



x)p 



+ 



I ilr dB [ (x)o[(x)+(IB! 1 (.v)n 2 (.v)+vx 



(75) 



We now proceed to the minimization of the increment of Trjc^cr] within the coherent state ansatz along 
the same lines as the previous subsection. First we optimize the noise db on the normalization factor II: 



db = -NU (^j dx dB x {x)(j)l(x) + dB* 2 {x)<j) 2 {x] 



(76) 



This choice leads to a vanishing noise term in Eq.(^). We insert this expression for db in the validity 
conditions Eq.([50|) and Eq.(^) and we obtain: 



Fa(x) = ~ 

in 



h + N I dx'V(x- x')\cp a (x')\ 2 



4>a(x) 



(77) 



Finally minimization of the contribution of the noise terms dB a with the constraint Eq. (|52|) is achieved 
with the choice 



dk 



f dt \ 1//2 

dBa{x) = {jn) + M J 



(V(k)) 



1/2 



(78) 



where the phases 9 a (k) are randomly generated as in Eq.(j7 ! 
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The first equation Eq.(|77|) fixes the deterministic evolution to the usual mean-field equation (|j). We 
note here that the mean- field term in Eq.(|77[) does not contain the normalization of the spatial density 
N\4> a (x')\ 2 by ||0 a || 2 , a feature present in the Fock state simulation (see Eq.(|20|)). This is a disadvantage of 
the coherent state simulation since this normalization factor appearing in the Fock state simulation has a 
regularizing effect: the norms \ \<pa\ \ may indeed deviate significantly from unity in the stochastic evolution. 
The second equation Eq.( |78|) determines the stochastic noise on the wave functions in a way very similar 
to the Fock state case Eq.(f38|). In particular the evolutions of </>i and 4>2 are still uncorrelated. The only 
difference is the disappearance of the projector Q a in the expression of the noise, which leads to an increased 
noise with respect to the simulation with Fock states. 

As in the previous subsection we now estimate the squared error A. We calculate the variation of Trfcr^a] 
for the choice of noise Eq.( 



1 dTr[a^a] NV(0) 



Trfcrta] dt 



a=l,2 



(79) 



The average over all stochastic realizations of the norm squared of the wave functions can be calculated 
exactly: 



i 2 ) 

/ stoch 



'nil / . (*) = e 



^(O)A/|| 0q ||2\ 



/ stoch ^ ^ 



This leads to a remarkable identity on the trace of a' a: 

lnTr[crtcr]\ (t) = (lnTrfcrtfjA (0) + N (e tVi ° )/h - l) ( V 

' stoch > ' ot""^ V / \ ^— 



stoch 



a=l,2 



^ ^ ^ stoch ^ ^ 



(80) 



(81) 



Using finally the concavity of the logarithmic function, leading to the logarithm of a mean being larger than 
the mean of the logarithm we obtain a lower bound on the squared error A on the N— body density matrix: 



A(t) + Tr[p 2 ] > .lexp \ 2BN [e tv ^l h - 1 
where we have introduced the constant quantities 

A = exp [/lnTr^cr]) _ (0) 



stoch 



B 



Q = l,2 



I 2 ) 



stoch 



(0) 



(82) 

(83) 
(84) 



It is quite remarkable that in the limit of times shorter than h/V(0) this lower bound scales exponentially 



with time as the upper bound Eq.(74) obtained for Fock states. Consequently the simulation scheme with 
Fock states is likely to be more efficient that the simulation with coherent states. This will indeed be the 
case in the numerical examples given in the next sections. 



3.3 The constant trace schemes 



We have given the expression of the one-body density matrix in terms of 4>a(x) for the simulation with 
Fock states Eq.([l3|). This expressions shows that is very sensitive in the large N limit to fluctuations of 
(02 l^i) - The same remark applies to two-body observables. In order to improve the statistical properties of 
the simulation one can consider the possibility of a simulation scheme with (4>2\<Pi) = 1 at any time. This 
actually corresponds to a conserved trace of each single dyadic a(t). This possibility is analyzed in § 3.3.1 ; it 
is extended to the coherent state simulation in 



3.3.2, leading to the well-known positive P-representation 



formalism. 
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3.3.1 Simulation with Fock states 



Within the Fock state ansatz, the conservation of the trace of the dyadic Tr[cr] can be achieved by (i) 
choosing the transverse noises dB^ according to the formula Eq.(|68[) and (ii) using the expression Eq.(|65|) 
for the longitudinal noise with X = 1. Point (ii) gives: 



d lx = -((^i)- 1 J dx 4>* 2 (x)dBi(x). (85) 
The forces terms X a and F a are fixed by the conditions (|39|)-([4~l]): 



N - 1 



Ai dt = -^Y-^\^y 2 I dx dx ' 4>l{x)4> ¥ 2 {x')dB^{x)dB^{x') (86) 

and 



F^(x) dt = (N- l)<0 2 |0i> _1 J dx' <j)* 2 (x')dB^(x')dB^-(x). (87) 
The expressions for dj2> ^2 and F 2 (x) are obtained by exchanging the indices 1 and 2 in these results. 



3.3.2 Simulation with coherent states 



In the case of the coherent state ansatz, the value of da s , which is the zero-mean noise term entering the 
variation of the dyadic a during a time step dt, is given in Eq.(^0|). The requirement of a constant trace 
Trfc] = rie^^ 2 '^ 1 ^ leads to the following condition on the noise terms 



db + NU / dx {(^{xjdB^x) + dBl{x)<j)x{x)) = 0. 



We choose the noise terms dB a as in Eq.(78). The remaining parameters F a are now unambiguously 
determined by (|50|)-(^T1): 



in 

F2(x) = \ 
in 



h + N j dx' ^(x^Vix-x'^x') 
h + N [ dx' (pl(x')V(x - x')4> 2 {x') 



4>i{x) 
<fo{x). 



(89) 
(90) 



This scheme exactly recovers the stochastic evolution in the positive P-representation, which was originally 
obtained with a different mathematical procedure p3f|. 



4 Stochastic vs. exact approach for a two-mode model 

In order to test the convergence of the stochastic schemes developed in the previous section we now apply 
this method to a simple two-mode system for which the exact solution of the iV-body Schrodinger equation 
can also be obtained by a direct numerical integration. This allows (i) to check that the stochastic methods 
when averaged over many realizations give the correct result indeed, and (ii) to determine the statistical 
error for each of the four implementations of the stochastic approach ('constant trace' vs. 'simple', Fock vs. 
coherent states). 

The toy-model that we consider describes the dynamics of two self-interacting condensates coherently 



coupled one to the other. It can be applied to the case of two condensates separated by a barrier [24 



(Josephson-type coupling) or condensates in two different internal states coupled by an electromagnetic 
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field |25|] (Rabi-type coupling). In this model we restrict the expansion of the atomic field operator to two 
orthogonal modes, 

ip(x) = au a {x) + bu b (x). (91) 
The Hamiltonian Eq.(|l]) takes the simple form: 

n = ^ ( 6 ts + Sta) + Hk (a) 2 a 2 + IH 2 ) (92) 

where a, b annihilate a particle in modes u a and Ub, K characterizes the strength of the atomic interactions 
inside each condensate and f2 is the Rabi coupling amplitude between the two condensates. Here we have 
restricted for simplicity to the case where (i) the condensates have identical interaction properties, (ii) the 
interactions between atoms in different wells are negligible, and (iii) the Rabi coupling is resonant. The 
most general two- mode case could be treated along the same lines. 

The direct numerical solution of the Schrodinger equation is performed in a basis of Fock states \n a ,nb) 
with n a f, particles in modes u a ,Ub- The numerical integration is simplified by the fact that n a + n b is a 
quantity conserved by the Hamiltonian evolution. We start with a state in which all atoms are in mode 
Ub, either in a Fock state \n a = 0,nb = N) (for the Fock state simulations) or in a coherent state oc 
exp(A^ 1 / 2 6^)|0, 0) (for the coherent state simulations). We watch the time evolution of the mean fraction of 
particles in mode u a , p a = (a'a)/N. 

Mean-field theory (the Gross-Pitaevskii equation), valid in the limit N 3> 1 with a fixed kN/Q p6[ , 
predicts periodic oscillations of (a^a)/N; the peak-to-peak amplitude of the oscillations is equal to unity 
if kN/Q < 1, and is smaller than one otherwise 27]. In the exact solution the oscillations are no longer 



periodic due to emergence of incommensurable frequencies in the spectrum of H. 

In the simulation method we evolve sets of two complex numbers representing the amplitudes of the 
functions (j>i(x) and 4>2(x) on the modes u ay b{x) (plus the n coefficient in the coherent state case). The 
results are presented in figure [l] for N = 17 particles and n/Q = 0.1, together with the result of the direct 
integration of the Schrodinger equation. 

The first row in the figure concerns the constant trace simulations. Figure |]a shows results of the 
simulation based on the positive P-representation, that is the constant trace simulation with coherent 



states. As well known [17] this scheme leads to divergences of the norm \\4>i\ \ \ \4>2\\ for some realizations of 
the simulation. We have cut the plot in figure |] at the first divergence. The same type of divergences occurs 
in the constant trace simulation with Fock states (figure |l]b). Note however that the characteristic time 
for the first divergence to occur is somewhat longer. We have checked for these constant trace simulations 
that the probability distribution of \\4>i\ \ WfoW broadens with time, eventually getting a power law tail. The 
corresponding exponent a decreases in time below the critical value a cr it = 3 for which the variance of 
ll^i]] 1 1 <p2 1 1 becomes infinite. This scenario is identical to the one found with the positive P-representation 
0- 

The simple simulation schemes plotted on the second row of figure [l| provide results which are at all time 
in agreement with the direct integration within the error bars. Contrarily to the constant trace schemes we 
do not observe finite time divergences in the simple schemes. For a given evolution time we have checked 
that the error bars scale as 1/a/AA where Af is the number of stochastic realizations. For a given Af we found 
that the error bars increase quasi-exponentially with time. 

The noise in the simple simulation schemes is investigated in more details in figure ^ which shows the 

error estimator /Trffj^cr] ) as function of time, for coherent states in figure |2]a and for Fock states in 
n stoch 

figure |2p. The coherent state result confirms the prediction Eq.(p2|). The Fock state result is found to 
be notably smaller than the upper bound Eq. (|74}) . This is due to the fact that the terms proportional 
to (<p a , <fioi\V\<fiai <Pa) i n Eq.(|73]) are not negligible as compared to the term V(0). We have checked these 
conclusions for various values of N and n/Sl. 
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Figure 1: In the two- mode model mean fraction of atoms in the mode u a as function of time, obtained with (a) the 
positive P-representation, (b) the Fock state simulation with constant trace, (c) the simple simulation with coherent 
states, and (d) the simple simulation with Fock states. The solid line represents the average over TV = 2 x 10 5 
simulations, with corresponding error bars. The dashed line is the direct numerical solution of the Schrodinger 
equation. The number of atoms is N = 17, initially all in mode u^. The interaction constant is k — 0.117. The time 
step used in the numerical stochastic calculation is £ldt = 10~ 3 . The calculations in (a) and (b) have been stopped 
after the divergence of one realization. 




Figure 2: Statistical error on the TV— body density matrix for the two-mode model: (a) 'simple' scheme with coherent 
states and (b) 'simple' scheme with Fock states. The solid line is the numerical result of the simulations. The dashed 
lines in (a) and (b) correspond respectively to the lower and upper bounds Eq.(|82]) and Eq.(|74|). The parameters are 
the same as in figure [|. 
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For a large number of particles it is known |2S] that the oscillations of (a^a) experience a collapse 
followed by revivals. These revivals are purely quantum phenomena for the field dynamics and they cannot 
be obtained in classical field approximation such as the Gross-Pitaevskii equation. We expect to see a 
precursor of this phenomenon even for the small number of particles N = 17. As the simple scheme 
simulation with Fock states is the most efficient of the four schemes for the investigation of the long time 
limit, we have pushed it to the time at which a "revival" can be seen, as shown in figure || This figure is 
obtained with M = 10 8 simulations. 



0.1 



0.05 




Figure 3: Fraction of atoms in mode u a in the two-mode model, for the parameters of figure [j]. The dashed line is 
the direct numerical solution of the Schrodinger equation. The solid line with error bars is the result of the 'simple' 
scheme simulation with Fock states with TV = 10 8 realizations. To keep a reasonable computation time with such a 
large value of AT we have increased the time step in the numerical stochastic integration by a factor 25 with respect to 
figure [j]. This explains the small systematic deviation of the simulation result from the exact result visible for example 
at time fit = 2.6. The quantum phenomenon of collapse and revival of the oscillation amplitude clearly apparent on 
the exact result is well reproduced by the simulation. 



5 Stochastic approach for a one-dimensional Bose gas 



The interacting Bose gas is in general a multi-mode problem, and the simulation schemes may have in this 
case a behavior different from the one in a few-mode model such as in §||. We have therefore investigated 
a model for a one-dimensional Bose gas. The atoms are confined in a harmonic trap with an oscillation 
frequency to. They experience binary interactions with a Gaussian interaction potential of strength g and 
range b: 



V(x - y) 



!J 



(27r)V2b 



exp [-( x - y f/(2b 2 )]. 



(93) 
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At time t = all the atoms are in the same normalized state <j) solution of the time independent Gross- 
Pitaevskii equation 



H4>{x) 



2m dx 2 



+ -mco 2 x 2 c 



x) + (N — 1) / dx' V(x - x')\0(x')\ 2 (j)(x) 



(94) 



At time t = + the trap frequency is suddenly increased by a factor two, which induces a breathing of the 

cloud simp. 



This expected breathing is well reproduced by the numerical simulations. The mean squared spatial 
width R 2 of the cloud as function of time is obtained by taking O = Y^k=i ^i/N in Eq.(ll) where xt is the 



position operator of the k— th particle. The quantity R 2 is shown in figure Q for the simulation schemes with 
Fock states. One recovers the key feature of the constant trace simulation, that is a divergence of the norm 
\\4>i\ \ 1 1 ^2 1 1 in finite time for some realizations. Before the occurrence of the first divergence the stochastic 
variance of the size squared of the cloud, defined as 



Var(P 2 ) stoch 



1 

Jf 



£ [R 2 (t) - R 2 (t)]' 2 with Rl(t) =Ke (N: <${t)\6\N : ^\t)) 



(95) 



is notably smaller in the constant trace scheme than in the simple scheme, as shown in figure [|a. This 
contrast between the two schemes for the statistical error on one-body observables was absent in the two- 
mode model of §£|. 

We have also investigated the noise on the iV— body density matrix characterized by ( Trfcr^o"] ) (see 

\ / stoch 

figure ||b). As expected this error indicator is smaller with the simple scheme. For this simple scheme it 
varies quasi exponentially with time with an exponent 7 ~ 4w, which is smaller by a factor roughly 2 than 
the one of the upper bound Eq.(74). This difference is due to the fact that the range b of the interaction 
potential is chosen here of the same order as the size R of the cloud so that the terms (cf) a , (f> a \V\(f> a , a ) 
neglected in the derivation of the upper bound are actually significant. We have checked for various ranges 
b much smaller than R that 7 then approaches the upper bound 2NV(0)/h. 



6 Conclusion and perspectives 

In this paper we have investigated a general method to solve exactly the N— body problem in the bosonic 
case. The principle of the approach is to add to the usual mean-field Gross-Pitaevskii equation a fluctuating 
term. We have determined the general conditions ensuring that the average over all possible realizations of 
this stochastic equation reproduces the exact N— body Schrodinger equation. 

This idea already received a particular implementation in quantum optics, in the frame of the positive 
P-representation. We recover here the scheme based on the positive P-representation as a particular case 
of a simulation evolving coherent states of the bosonic field with the constraint that the trace of the density 
operator should remain exactly equal to unity for each single realization. This provides a simple derivation 
of the stochastic evolution within this representation alternative to the usual one [^J based on analyticity 
properties. 

Among the many possible implementations of the general stochastic approach we have also investigated 
schemes evolving Fock states (that is number states) rather than coherent states. This is well suited to 
situations where the total number of particles is conserved. In particular we have identified a scheme 
preserving exactly the trace of the density operator which is for number states the counterpart of the one 
based on the positive P-representation. 

Schemes with constant trace are subject to divergence of the norm of some realizations in finite time. 



This effect, already known in the context of the positive P-representation 17], makes these schemes difficult 
to use. 
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Figure 4: Mean squared spatial width R 2 of a harmonically confined cloud of N = 10 atoms as function of time. The 
breathing of the cloud is induced by an abrupt change of the trap frequency from uj to 2u> . The width R is measured 
in units of the harmonic oscillator length aho = (S/fraw)) 1 ' 2 . The interaction potential is chosen such that b = 0.5ah o 
and g = OAhcuciho leading to a chemical potential \i = I.ITiuj in the Gross-Pitaevskii equation Eq. (|94|) . The calculation 
is performed on a spatial grid with 32 points ranging from — 6a^ to +6ay lo (with periodic boundary conditions). •: 
constant trace simulation with J\f = 1000 realizations. For ujt > 3.5 a divergence has occured for one of the realizations 
and the calculation has been stopped. □: simple scheme simulation with AT — 40000 realizations. 




Figure 5: For the one-dimensional Bose gas in the conditions of figure [|, (a) stochastic variance of the size squared 
of the cloud and (b) noise on the N— body density matrix. Solid lines: simple scheme with Fock states. Dashed lines: 
constant trace scheme with Fock states. 
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In order to overcome this divergence problem we have investigated schemes in which the condition on 
the trace is relaxed. We have chosen instead to minimize the statistical spread on the N— body density 
matrix, which gave rise to the 'simple' schemes, either with coherent states or Fock states. In this case the 
N— body density operator is obtained as a stochastic average of dyadics such as |coh : iV 1 / 2 0i)(coh : iV 1 / 2 ^! 
or \N : (f>i)(N : 02 1, where the evolutions of 0i and 02 are fully decoupled. The deterministic parts are given 
by Gross-Pitaevskii equations, which preserves the norm of 01,2, contrarily to the case of constant trace 
schemes. The decoupling between the evolutions of 0i and 02 allows a reinterpretation of our representation 
of the N— body density operator. If the initial density operator is given by p(t = 0) = \N : 4>q)(N : 0o| it 
will be given at time t by 



with the iV— particle state vector 



p(t) = ]*(*)> W)l (96) 



1 M 

\#(t))=]hn^ O0 jf y £ t \N:<l>V\t)). (97) 

In this expression (jyi) are stochastic realizations with the initial condition (fyj)(t = 0) = 0o- 

The 'simple' schemes have much better stability properties than the constant trace schemes: differ- 
ently from the case of constant trace schemes, the deterministic evolution of the 'simple' schemes has a 
Gross-Pitaevski form and thus conserves the norms ||0i,2||- This condition, together with the upper bound 
Cl,2 < V(0)dt\ |0i,2| | 2 /ft on the eigenvalues Ci,2 of the noise covariance operator dB a (x)dB^(x') , can be used 
to prove that the stochastic equations possess a finite, non-exploding solution valid for all times (see [19] in 
page 94, ||| in §4.5). 

We have numerically applied the simulation schemes to a two-mode model and to a one-dimensional 
Bose gas. In both cases we found that the constant trace schemes lead to some diverging realizations, while 
the simple schemes lead to a statistical spread on the N— body density operator increasing exponentially 
with time with an exponent 7 oc NV(0). The simple schemes are therefore not well suited to determine 
small deviations from the mean-field approximation in the large N limit but can be more efficiently applied 
to systems with a small number of particles, such as small atomic clouds tightly trapped at the nodes or 
antinodes of an optical lattice. 

In the one-dimensional numerical example of this paper we have presented results for a simple one-body 
observable, the size of the atomic cloud. We have actually extended the calculations to more elaborate 
observables such as the first order and the second order correlation functions of the field. We have not 
presented the results here as the initial state of the gas was taken to be a (not very physical) Hartree-Fock 
state. We are presently working on the possibility to generate a more realistic initial state such as a thermal 
equilibrium for the gas, by extending our stochastic approach to evolution in imaginary time. 

This work has several possible perspectives of extension. One can first use as building block a more 
sophisticated ansatz than the Hartree-Fock state \N : 0), such as Bogoliubov vacua (that is squeezed states 
of the atomic field) or a multimode Hartree-Fock ansatz (that is an arbitrary coherent superposition of 
number states in several adjustable modes of the field). One can also look for approximate rather than 
exact stochastic solutions to the N— body problem but that would be better than mean-field approaches in 
some given situations. We hope to address some of these perspectives in the near future. 
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